Nonlinear quantum heat transfer in hybrid structures: Sufficient conditions for 

thermal rectification 
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We present a unified description of heat flow in two-terminal hybrid quantum systems. Us- 
ing simple models, we analytically study nonlinear aspects of heat transfer between various reser- 
voirs: metals, solids, and spin baths, mediated by the excitation/relaxation of a central (subsystem) 
mode. We demonstrate rich nonlinear current-temperature characteristics, originating from either 
the molecular anharmonicity, or the reservoirs (complex) energy spectra. In particular, we establish 
sufBcient conditions for thermal rectiflcation in two-terminal junctions. We identify two classes of 
rectiflers. In type-A rectifiers the density of states of the reservoirs are dissimilar. In type-B recti- 
fiers the baths are identical, but include particles whose statistics differ from that of the subsystem, 
to which they asymmetrically couple. Nonlinear heat flow, and speciflcally thermal rectiflcation, 
are thus ubiquitous effects that could be observed in a variety of systems, phononic, electronic, and 
photonic. 

PACS numbers: 63.22.-m, 44.10.-|-i, 05.60.-k, 66.70.-f 



I. INTRODUCTION 

Understanding heat transfer in two-terminal hybrid 
structures is of fundamental and practical importance, 
for controlling transport at the nanoscale, and for realiz- 
ing functional devices [l|, Among the systems that 
fall into this category are metal-molecule-metal junc- 
tions, the basic component of molecular electronic de- 
vices. Here the excessive energy generated at the molec- 
ular core should be effectively removed for realizing a sta- 
ble device, as demonstrated theoretically [3] and exper- 
imentally 0, [H| • Another composite structure of funda- 
mental interest is a dielectric-molecule-dielectric system, 
where vibrational energy flowing between the compo- 
nents activates reactivity and controls dynamics 0, ■ 
Phononic junctions are also captivating for understand- 
ing the validity of the Fourier's law of thermal conduc- 
tion at the nanoscale . Single- mode radiative heat 
conduction between ohmic metals was recently detected, 
manifesting that photon-mediated thermal conductance 
is quantized [ll[. Other hybrid systems with interest- 
ing thermal properties are electronic spin- nuclear spin in- 
terfaces [13], metal-molecule-dielectric contacts [Ts'l and 
metal-superconductor junctions [14]. 

From the theoretical point of view, the basic challenge 
here is to understand the role of nonlinear interactions in 
determining transport mechanisms, and specifically, in 
inducing normal (Fourier) heat conduction, either within 
classical laws lol. [l5l. [l6l o r based on quantum mechanical 
principles [iSjIlllliniO, [II, ^] . Treatments employed 
vary. Ballistic heat transfer in phononic systems was ex- 
plored using the Landauer-scattering formalism and 
within the generalized Langevin equation [16, 24, 25]. 
For interacting systems there are only few analytical re- 
sults, and the techniques employed include the Kinetic- 
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Boltzmann theory [261 ] . mode coupling theory [2] 
the non-equilibrium Green's function meth od j2ll .l2 
classical [15| and mixed classical-quantum [311] molecular 
dynamics simulations, and exact quantum simulations on 
simplified models [3^ . 

The master equation formalism [33| is another useful 
tool for studying quantum heat transfer across nano- 
j unctions [13, [s^, [S]- We have recently adopted this 
method in the weak system-bath coupling limit, and ex- 
plored various aspects of heat flow in simplified toy mod- 
els 0, [13, [H, [1§] ■ Basically, we ask ourselves the fol- 
lowing question: What is the connection between the 
microscopic Hamiltonian and the heat current across the 
system? More specifically, how can we control the onset 
of nonlinear current-temperature behavior (stI . |38| , and 
what conditions should the Hamiltonian satisfy for the 
system to manifest normal conduction j20|? While our 
treatment has been typically limited to a minimal sub- 
system, the formalism has two main advantages: (i) Both 
harmonic and anharmonic systems can be treated within 
the same footings, (ii) Analytical results obtained can 
pinpoint on the underlying dynamics, typically obscured 
in numerical simulations. 

In this work we extend our recent letters [1^, [40| . 
and develop a unified description of temperature-driven 
quantum heat transfer at the nanoscale. We focus on 
two-terminal devices, including a central quantized unit 
(subsystem) and two bulk objects (referred to as termi- 
nals/contacts/reservoirs/baths), whose temperatures are 
externally controlled. Our description can be applied on 
several systems, including electron-hole pair excitations, 
phonons, photons or spins. However, our formalism does 
not describe heat flow due to the transport of hot charge 
carriers [1]. 

We analytically explore current-temperature charac- 
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teristics in various systems, and seek to connect the non- 
linear behavior to the microscopic parameters. We also 
investigate the temperature dependence of the thermal 
conductance, and observe rich behavior, depending of 
the details of the model. In particular, we discuss a sim- 
ple nonlinear effect, thermal rectification, an asymmetry 
of the heat current for forward and reversed tempera- 
ture gradients. This phenomenon has recently attracted 
considerable theoretical 0, [4iy4iE H, S [H, 113, 
m, m, [E^l and experimental[5iniiattention. While 
most theoretical studies have analyzed this phenomenon 
in phononic systems using classical molecular dynamics 
simulations, our prototype model can describe thermal 
rectification of several subsystems (vibrational or radia- 
tion modes) and reservoirs (spin, metal, dielectrics) at 
the same footing, see Fig. [T] Moreover, using our simple 
model we establish sufficient conditions for rectification 
[40j . We identify two classes of rectifiers: (i) "Type A", 
where the reservoirs are dissimilar i.e. of different den- 
sity of states (DOS), (ii) "Type B", where the baths are 
identical, but their statistics differ from that of the sub- 
system, combined with unequal coupling strengths at the 
two ends, as explained below. 

The content of the paper is as follows. In Section II we 
present our model. In Section III we derive equations of 
motion for the subsystem population using the quantum 
master equation approach. We consider several models 
for the subsystem and for the contacts, and derive ana- 
lytical expressions for the heat current through various 
conducting junctions. In Section IV we trace these junc- 
tions' (nonlinear) current-temperature characteristics to 
the Hamiltonian parameters. In particular. Section V is 
focused on a specific nonlinear effect, thermal rectifica- 
tion. We analytically identify two types of rectifiers, and 
demonstrate with numerical simulations the tunability of 
the effect. Conclusions follow in Section VI. 



operator J, 

J = Tr[Jp]; J^'-[Vl,Hs]+'-[Hs,VrI (3) 

where the expectation value is defined positive when the 
current is flowing L to R. The subsystem Hamiltonian 
assumes a diagonal form, and we also consider separable 
couplings 

Hs - ^^„|7i)(n|; 

71 

K = X.SB,; 5 = ^5„,,„|m)(7i|. (4) 

Here S" is a subsystem operator and B^, is an operator 
in terms of the bath degrees of freedom. For simplic- 
ity we set Sm.n = Sn,m- In what follows we consider 
situations where Bl and Br have the same dependence 
on the (respective) bath degrees of freedom, with differ- 
ent prefactors Xl ^ A^. We refer to this scenario as 
" parametric asymmetry" , rather than " functional asym- 
metry", resulting from dissimilar B's. Note that if the 
commutator [_ff5,S'] = 0, the heat current trivially van- 
ishes. 




FIG. 1: Examples of two hybrid systems treated in this work, 
(a) Single mode heat transfer between a solid and a spin bath. 
The central unit can represent either a vibrational or a radi- 
ation mode, (b) Phonon to exciton energy exchange. 



II. MODEL 



We consider generic hybrid structures where a central 
unit Hs interacts with two reservoirs (v — L, R) oi 
temperatures = (3~^ [kB = 1) via the coupling terms 

H = HI + HI + Hs + Vl + Vr. (1) 

The heat current from the u bath into the subsystem 
is given by ^, J, = §Tr ([ijo - i/g, K]p) ; {h = 1), 
where p is the total density matrix, and we trace over 
both the subsystem and reservoirs degrees of freedom. 
In steady-state the expectation value of the interaction 
is zero, Tr (^j^p) = 0, and we obtain 

JL^tTT{[VL,Hs]p): JR^tTv{[VB,Hs]p). (2) 

Since in steady-state the currents are equal, Jl = —Jr, 
we introduce a symmetric definition for the heat current 



III. DYNAMICS 

Employing the Liouville equation in the interaction 
picture, the elements of the total density matrix satisfy 

^ = ~^[V{tlp{0)U^ 

- / dr[y(i),[y(T),p(T)]]„,„, (5) 

Jo 

where V = Vl + Vr, and V{t) are interaction picture 
operators. Following standard weak coupling schemes 
|33| . we proceed by making four assumptions: (i) We 
first make the presumption that Tr[y(i), /9(0)] = 0, i.e. 
the mean value of the interaction Hamiltonian, averaged 
over the initial density matrix, is zero, (ii) We factorize 
the density matrix of the whole system, at all times, by 



3 



the product p{t) = cr(t)pL{TL)pR{T!i)- Here is inde- 
pendent of time, PuiT^) = e-^°^'^"/Tr^[e-""/'^"], and 
a{t) is the subsystem density matrix obtained by trac- 
ing out the reservoirs degrees of freedom from the total 
density matrix, a{t) — TrB[p{t)]. Trs denotes trace over 
both the L and the _R-baths degrees of freedom, (iii) 
We take the Markovian Umit, assuming that the reser- 
voirs' characteristic timescales are shorter tlian the sub- 
system relaxation time, (iv) We neglect contributions 
from quantum coherences, i.e. in the long time limit we 
assume that the nondiagonal elements of the system den- 
sity matrix vanish. This can be justified in the weak cou- 
pling limit adopting the initial condition fTn^m(O) ~ 0. 
Following these assumptions, the Pauli master equation 
for the populations P„(i) = TrB[pn,n(i)] is obtained [ssj . 

- Pnit)J2\S„.,n\'K^raiT.)- (6) 

The Fermi golden rule transition rates are given by 



Here E„ 
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(7) 
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T^"c,y[p^(Ti,)Bi,{T)BT-(0)] is the correlation function of the 
h' environment. In steady-state P„ = 0, and we normalize 
the population to unity P„ = 1. It is straightforward 
to show that under ^ the steady-state current ^ be- 



n,m 

(8) 

with the population obtained by solving ^ in the long 
time limit. For details see Appendix A. Our description 
to this point is general, as we have not yet specified nei- 
ther the subsystem nor the terminals. 



A. Specific models for the subsystem Hamiltonian 

We consider two representative models for the subsys- 
tem Hamiltonian and its interaction with the baths. In 
the first model the subsystem is a harmonic oscillator 
(HO) of frequency oj, Hs = nLo\n){n\ This can de- 
scribe either a local radiation mode 'SSl , or an active 
vibrational mode of the trapped molecule [38] . For the in- 
teraction operator we assume S — ^/n\n) (n — 1 1 -I- c.c, 
motivated by the bilinear form oc xi?^, a; is a subsys- 
tem coordinate [s^. This implies that only transitions 
between nearest states are allowed, 



J — oo 

K^i^niT.) = e-P--K^^_,{T,). (9) 



We also introduce the short notation 



(10) 



where f^, defined through encompasses the effect of 
the bath operators. Solving ^ in steady-state using the 
above expressions for the rates, the heat current ^ can 
be analytically calculated [sj, 



j(HO) ^ 



n%{--Lo)/kL{TL)+n'^{-Lo)/kR{TRy 



(11) 



where ng{u;) — ^e^"'^ — l] is the Bose-Einstcin distri- 
bution function at T^, — 1/(3^. 

Our second subsystem is a two-level (spin) system 
(TLS). Here Hs = ^cr^, and we employ a nondiagonal 
interaction form S = Ux- These terms can represent an 
electronic spin rotated by the environment 12]. It can 
also describe an anharmonic (truncated) molecular vibra- 
tion that is dominating heat flow through the junction 
[s^,!!^]- Re-calculating the long-time population ([6]), the 
heat flux reduces to 



j{TLS) 



Cj[7lf(w) -7lf(w)] 



n^(-c.)/fc^(rL) + n'l{-Lo)/kR{TR) 



(12) 



with the rates ([9]) and the spin occupation factor ng{Lu) = 

Expressions pip and demonstrate that in the 

weak coupling limit the effect of the environment enters 
only through the relaxation rates k"^ , evaluated at the 
subsystem energy spacing lo. In Section IV we further 
introduce a three-level system, an intermediate structure 
between a two-level (strongly anharmonic) system and 
an harmonic object. More generally, given the subsys- 
tem anharmonic potential and an interaction operator 
S, one should proceed by (i) calculating (analytically or 
numerically) the vibrational spectra \n), (ii) evaluating 
the matrix elements of S, (iii) acquiring the steady-state 
levels populations by solving ([5]), and (iv) arriving at the 
heat current using 



B. Calculation of the rate constant 

We give next the explicit form for the relaxation 
rate for various physical thermal baths. The rate con- 
stants © induced by e.g. the L contact is given by 
the Fourier's transform of the bath correlation function 
(BL(r)Bi(0))j,^, with Bl = Y.Bli,\l){l'\ and Hi = 
^£^;]/)(Z]. Here \l) are the many-body states of the L 
reservoir with energies Ei . The rate ([9]) can be integrated 
to yield 



k'^{T,)^2^\lY,\Blu'\ 5{Eu-Ek 

k,k' 



- LO 



ZAP.) 



(13) 
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with Z^{j3^) = e~^"^'' , the partition function of the 
V bath; k, k' € v. 

(i) Distinguishable noninteracting particles. This en- 
vironment includes a set of independent (p — 1,2, 
particles. The Hamiltonian is given by summing all sep- 
arate contributions 

p p 

Within these terms, the relaxation rate reduces to 

P i,j 

g-/3..ep(j) 

X (^(£p(0-ep(j)+^) z {P ) ' (^^^ 

with the p-particle eigenstates \i)p and eigenvalues ep(i). 

Zp{(3w) = X^i e"''"'''^*' is the p-particle partition function. 
Specifically, for a bath of noninteracting spins we get 

r(r.) = r^(c.K(-^), (16) 

with the spin occupation factor n'g{uj) — [e^"'^ + l] ^ 
and the temperature independent coefficient 

r^H = 2nXl I (0|p \^6{lu + 6p(0) - £^(1)). 

p 

(17) 

(ii) Solid/Radiation field (harmonic bath). This bath 
includes a set of independent harmonic oscillators, cre- 
ation operator j . System-bath interactions are further 
assumed to be bilinear, 

i/0 = ^^,a^a,,,; = ^(a.^ + a^). (18) 
j j 

This leads to the relaxation rate 

r(r.) = -r^(^)n^(-c.), (i9) 

with the Bose-Einstein function n^(ti>) = [e^"^ — 1]^-'^ 
and an effective systcm-bath coupling factor 

T^L^) - 27rXlY,S{uJj-uj). (20) 

3 

(iii) Fermionic bath: metal. Consider a metallic terminal 
including a set of noninteracting spinless electrons, cre- 
ation operator j , where only scattering between elec- 
tronic states within the same lead are allowed, 

^° = E^4.c-^- s^ = E<«^-.^- (21) 

The transition rate between subsystem states (|T3]) can be 
written as [l^ 

k-'in) - -27rAX(-w)E'5(^'-^J-+^) 

hi 

xK(e,)-"F(e»+^)], (22) 



with the Fermi-Dirac distribution function np{e) = 
[g/3^(e-A'v) _|_ 2^]-!^ jg |;i^g chemical potential. One could 
also write 

r(r,) = -A"^(T,,u;)n^(-tu), (23) 

where 

A'in^iu) ^2Tr J de K(e) - n^(e + w)] F,(e). (24) 

The function F^{e) = ^(5(6 - -I- w)(5(e, - e) de- 
pends on the system-bath coupling elements and the spe- 
cific band structure. Assuming that the density of states 
slowly varies in the (subsystem) energy window the 
interaction function can be expanded around the chem- 
ical potential JJ.], F^{e) ~ F^{^^) + Iv^-^-^jj;^-, with 
a dimensionless number of order unity. Using this form, 
the integration in Eq. (|24p can be performed when the 
Fermi energy is much bigger than the conduction band 
edge ^i, ^ Ec- One can then write 

A''(r,,c^)«r^(c^) (^l + ,5,^^ , (25) 

where S is a constant of order one, measuring the devia- 
tion from a flat band structure near the chemical poten- 
tial [1^. The coupling constant is given by 

r^(Lj) = 27rtjF,(^,). (26) 

When = we find that A''{Ti,,uj) = r^(w), a temper- 
ature independent constant, and the harmonic limit (jl9p 
is retrieved [Hj. 

To conclude this discussion, assuming different types 
of reservoirs and system-bath couplings, we derived here 
three expressions for the relaxation rate ([T5)) [or equiva- 
lently dH)]: Eq. ([16]) assuming a spin bath, Eq. (flQ)) for a 
phononic environment with a linear coupling to the sub- 
system, and Eq. ([23| for a metallic bulk with electron- 
hole pair excitations 

{n'c;{~uj)Tg{u!); noninteracting spins 
— n^(— w)r^(a;); phonons; linear coupling 
-n^(-w)(l + (5^^)r^(w) Metal 

(27) 

In each case the coefficient F"^ reflects the system-bath 
interaction strength, whereas the temperature dependent 
function encloses the reservoirs properties. 

IV. NONLINEAR CURRENT-TEMPERATURE 
CHARACTERISTICS 

Based on the two models for the subsystem (HO, TLS) 
and the different types of baths, we can construct sev- 
eral two-terminal heat-conducting junctions. We con- 
sider next few examples, and manifest nonlinear current- 
temperature behavior. We will also introduce a three- 
level system (3LS), an intermediate structure between a 
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two-level (strongly anharmonic) system and an harmonic 
object. Formally, if we write the heat current as 



J{Ta,AT) = Y,(^kiTa)AT'' 



(28) 



with AT = Tl — Tr; Ta = T^+Tji, we ask ourselves what 
are the microscopic parameters that are incorporated in 
the nonlinear coefficients ak', k > 1, and how can we 
control the magnitude of these terms. We also define 
the thermal conductance as /C = limAT^o J / AT , and 
examine its temperature dependence. 

1. Harmonic Junction. Our first example is a fully 
harmonic system, incorporating two harmonic reser- 
voirs connected by a harmonic link, modeling vibra- 
tional/photonic heat transfer between two solids/ohmic 
metals [111, 113: assuming that a specific harmonic 
mode of frequency uj dominates the dynamics. The cur- 
rent is given by (fTTjl with the rates (fT9|) , resulting in the 
thermal Landauer expression [23l |. 



r r 

Here Tb — ri^^rn ^ (temperature independent) trans- 
mission coefficient, with the coupling elements T^ cal- 
culated at the molecular frequency lo. This equation 
describes ballistic thermal transport, where energy loss 
takes place only at the contacts. Note that in the clas- 
sical limit J oc AT, K. = Tb, and the current does not 
directly depend on the molecular (subsystem) vibrational 
frequency. For a harmonic junction nonlinear effects 
are thus purely quantum, originating from the quantum 
statistics. 

2. Spin Junction. Our second example is a spin-TLS- 
spin system representing e.g. a central electron interact- 
ing with two nuclear-spin environments under an applied 
magnetic field [12]. The current is given by Eq. p2| with 
the rates (|16|) . bringing in a spin-Landauer formula. 



J = Tblj [n|(cj) - ng{uj)\ 



AT-— AT3 + 0(AT5) 



(29a) 
(29b) 



J = Tsuj [ns{uj) 

^ a 



E 



n=l,3,5 



AT 

tT 



(30a) 
(30b) 



r r 

with the transmission coefficient Ts — p/, , To is eval- 

uated at uj, the central spin energy spacing. Unlike (j29bp . 
the linear response term here does depend on the sub- 
system frequency, /C = Tslu'^/T^, and nonlinear terms 
persist even in the high temperature limit. Note that 
both Eqs. (|29ap and (|30aP are in the form of Landauer 
formula, since only elastic scattering processes are effec- 
tive (system and reservoirs are identical). Our next two 
examples bring in deviations from the Landauer picture. 

3. Harmonic baths-spin subsystem junction. Consider 
again vibrational heat flow. However, unlike the fully 
harmonic model resulting in (I29p . the central unit is as- 
sumed here to incorporate anharmonic terms [stI . [39l . [50| . 



This can be implemented by modeling the subsystem by 
a truncated harmonic oscillator, for example a spin. The 
current across the device is evaluated using (fT2|l with the 
rates (flT 



J : 



LO 



r|ri(n|(^) 



/AT 



w/3„<l ^ \ - 

V T 

n=l \ 



-rg(l + 2ng(c.)) 



(31a) 
(31b) 



where xb = pf , pg measures the spatial asymmetry. 

Note that we could still organize Eq. piap in the form 
of the Landauer formula with a temperature dependent 
transmission coefficient [2l|. The high-temperature lin- 
ear response limit yields JC = TsLo/Ta. Nonlinear terms 
survive in (j31bp only due to the asymmetric coupling, 
since for xb — 0, J AT. Thus, quite interestingly, this 
anharmonic junction conducts linearly (in the classical 
limit), as long as it is fully symmetric. 

In Fig. [2] we exemplify the properties of the junctions 
([29|) - ([3T|) with asymmetric couplings F-'^ F^, taking uj, 
the subsystem characteristic energy, to be either of the 
order of the reservoirs temperature, a; ~ T^, or signifi- 
cantly lower. In the first case subplots (a)-(c) reveal that 
the current saturates at large AT, reflecting the quan- 
tum statistics. In contrast, in the high temperature limit, 
while a pure harmonic system (d) shows a linear current- 
AT characteristic, the other two systems (e)-(f), incor- 
porating some nonlinearities, reveal nonlinear behavior. 



(a) 






HO junction 





1 2 


(d) 






''ho junction 



^ 0.2 



10 




\-Tr 



FIG. 2: Nonlinear heat flow in various hybrid junctions, (a)- 
(c): uj = l,Ta=2. (d)-(f): uj = 0.1, Ta = 15. (a) and (d) are 
harmonic junctions with the current (|29[1 : (b) and (e) are spin 
junctions obeying (|30|) : (c) and (f) are harmonic baths-spin 
junctions obeying ((3l|- F^^IO, = 1 in aU flgures. The x 
axis is AT. 

4. Three-level junctions. In order to further eluci- 
date the role of a uniform energy spectra in nonlinear 
transport, we construct next a three-level system (3LS), 
an intermediate structure between a two-level (strongly 



6 



anharmonic) subsystem and an harmonic object. We as- 
sume that both the subsystem and the reservoirs include 
(equivalent) 3LS particles, Hs = J2n=i 2 3 En\n){n\; S = 
|1)(2| + |2)(3| + h.c. The current ^ reduces to 



), (32) 



^2 

with uji = E2 — El] LU2 = E3 ^ E2, and 

Di = A:i^2fe^3/fc3^2 + ^1^2 + ^2^1; 

D2 = k3^2k2^l/kl^2 + k3^2 + ^2^3- (33) 

Here ki^j — kf'_^j + k^-, with the rates defined in ([7]). 
Taking u = uji ^ lu2, and assuming the reservoirs in- 
clude collections of 3LS noninteracting particles of equal 
spacing to, we obtain the following expression, using p5|) 
(c^ < T„), 



9 Ta 



Ti=l,3,5... 



(34) 



with Ts = and = 2nXlj:j^\{i\b,,p\j)p\H{Lu + 

tp{i) — £pij))', i < j- This expression is essentially similar 
to pop . It is notable that only odd terms survive, irre- 
spective of the spatial asymmetry (F^ ^ F^). It can be 
shown that even terms participate only when the energy 
spectra along the device becomes inhomogeneous. 

5. Energy flow between metals. Energy transfer be- 
tween metals, mediated by the excitation of a single ra- 
diation mode, has recently attracted considerable exper- 
imental and theoretical interest [111, [3^, [s^, H^l- Anal- 
ogous junctions are the core of molecular electronics, 
where the bridging modes are the vibrations of the 
trapped molecules. For example, heat dissipation in a 
metal surface-Cgo-STM junction was demonstrated to be 
controlled by vibrational decay into the tip, generating 
electron- hole pair excitations [1, [l^. Heat conduction 
through a DNA-gold composite, suspended between two 
electrodes, was found to be dominated by phonon trans- 
port m. 

When both metals are ohmic, the current in the weak 
coupling limit is given by ([29)1 [39|. As we show next, 
for structured reservoirs nonlinear effects emerge. We 
calculate the heat current using Eq. pT|) with the rates 
assuming /i = /i^ and 5 — 5„, i.e. the reservoirs 
are equivalent, and are maintained at the same chemical 
potential. This leads to 



J ^ bJ 



rf r|^(l + ^^)(l + 5^) [n'i.iio) ~ n^Bi^)] 



(F^ + Tfy 



(35a) 



(35b) 



In deriving (j35bp we assumed that STa/^^ < 1. Here 
Xf = pL I measures the asymmetry in the system- 

p-L pi? 

bath coupling and Tp = jTrrrFR- The following observa- 

tions can be made: (i) The nonlinear contributions are 
all related to a finite S, measuring the deviation from a 
constant density of states 39]. If (5 = 0, the harmonic 
limit is recovered, (ii) The existence of even terms, e.g. 
a AT^ term, requires some asymmetry, xf 0, as we 
discuss below, see Eq. pS]) . (in) The conductance of 
the junction scales like K, — Tp(l -I- STa/Iji), in sharp 
contrast to the behavior of phononic systems. 

In Fig. [3] we study energy transfer between metals 
employing a Lorentzian density of states of width 7, 



Du{e) 



7/27r 

— w 



We observe a small 'negative differ- 



ential conductance', a decrease of the current with in- 
creasing temperature different at large temperature bias. 
The effect prevails when 7 < cj, i.e. the reservoirs density 
of states is notably changing over the subsystem energy 
scale. The data was generated by employing (fTT|) with 
the rate The inset presents the F function ^I^, 

evaluated with the Lorentzian DOS. 

Summarizing, while in pure harmonic systems non- 
linear dynamics is linked to quantum effects (|29p . spin 
systems manifest nonlinear dynamics irrespective of the 
subsystem frequency, ([30]) and (p4)l . In a simple model 
for anharmonic vibrational heat fiow (|3ip nonlinearities 
are attributed to the spatial asymmetry, while nonlin- 
ear effects in excitonic energy transfer are linked to the 
metals' energy dependent density of states ((35l) . The con- 
ductance temperature dependence also varies. For har- 
monic systems it is independent of temperature, while for 
anharmonic phononic systems it decays like 1/T, in gen- 
eral agreement with experimental results [56j . In metal- 
single mode-metal junctions the conductance increases 
with temperature, due to the enhancement of the transi- 
tion rate with T. 



V. THERMAL RECTIFICATION: SUFFICIENT 
CONDITIONS 

We focus next on a specific nonlinear effect, ther- 
mal rectification, asymmetry of the heat current for for- 
ward and reversed temperature difference, |J(-|-Ar)| ^ 
|J(— Ar)|. Formally, the onset of this phenomenon is 
identified by the existence of even terms, ak=2n 7^ 0; 
n = 1, 2..., in the expansion 

We discuss here sufficient conditions for thermal recti- 
fication by analyzing the currents in (jlip and . These 
expressions are naturally not the most general results for 
heat transfer across hybrid structures. Rather, they relay 
on few assumptions as explained in Section HI: The sub- 
system and reservoirs are weakly interacting, the baths 
are assumed to be held at thermal equilibrium, and the 
dynamics is markovian. Furthermore, specific forms for 
the subsystem are employed. However, these analytical 
forms, derived from a prototype model, can still illumi- 
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FIG. 3: Negative differential conductance in a metal-single 
mode- metal junction with 7 = 0.2 (solid line). rj;^=r^=l, 
Ta = 3, tj = 0.5. The dashed curve was generated using 7 = 5, 
manifesting behavior similar to the fully harmonic case (|29|) . 
Inset: The function F{e) [see (|24p ] vs. energy for 7 = 0.2 
(solid line); 7 = 5 (dashed line). 



nate on the basic mechanisms involved in thermal rec- 
tification. Having said so, we return to Eqs. (fTTj) and 
P^ . and analyze their structure for forward (T^ = Th\ 
Tr = Tc) and reversed (T^ = Tc; Tr = Th) tempera- 
ture gradients, Th > Tq- We note that in each expres- 
sion separately these currents deviate in magnitude if the 
denominators fulfill 



n"{~uj) n^i-uj) n^{-uj) n"{-uj) 



where n'^ is either the spin occupation factor or the Bose- 
Einstein distribution. Rearranging this expression and 
making use of (fTO|) we get 



(37) 



In what follows we identify two classes of rectifiers. In 
type-A rectifiers system-bath interactions are equivalent 
at both contacts, but the reservoirs have distinct prop- 
erties. In type-B rectifiers the reservoirs are equivalent, 
but the subsystem statistics is distinct from the baths, 
combined with some parametric asymmetry. 



A. Type-A thermal rectifier 

The rate constants (|13p can be expressed in terms of 
the reservoirs density of states. For example, the rate 



induced by the L contact is given by 



fc^(T) = 27rAi^|B,^;,f <5(£;, c.) 



11' 



27rAi 



(38) 



Here ZUP) - Ei 



J dee ^'^DL{e) denotes the 



partition function of the L bath with the density of 
states DL{e) = 'Ei&{e — Ei). The function gL{e,uj) — 
J2i \B^{e,Ei)\'^5{e — Ei + oj) characterizes system-bath 
interactions. Taking Xl = Xr, we note that since the 
two sides of (|37p depend on different temperatures, the 
system rectifies heat if 



(39) 



besides special points in the parameter space, depending 
on the details of the model. Using ([55]) . this condition 
reduces to 



J dee-^'WL{e)gL{e,Lu) ^ J dee-P'Dn{e)gR{e,Lo) 



Jdee-0^DL{e) 



Jdee-0WR{e) 



(40) 



We next further assume that gL{e,uj) = gR{e,uj), i.e. 
system-bath interaction matrix elements are equal at 
both ends. Hence, the inequality ([M)) is satisfied when 



DL{e) ^ Dnie). 



(41) 



We thus recover a sufficient condition for thermal rectifi- 
cation: The density of states of the reservoirs should be 
distinct. Note that at least one of the reservoirs should 
have an energy dependent DOS. If both reservoirs are 
harmonic, D^j — c^, with Ci, a constant, rectification is 
absent even for cl cr. It is thus sufficient to have one 
of the reservoirs incorporating anharmonic interactions. 
If both reservoirs are nonlinear, rectification takes place 
if the energy spectra are distinct. This discussion relays 
on the assumption that the spectral function gje, ui) ex- 
plicitly depends on energy, see Eq. ([iU)) and pi]. 

Going back to the condition (HO]), we Taylor-expand 
the interaction function around lu, gy{e,uj) ^ ai(cj) + 

Q!2(w)e, 



j dee-f^''DL{e)e J dee-'^'DR{e)e 
Jdee-P^DL[e) ^ J dee~P^DR{e) ' 



(42) 



We identify the left (right) hand side as the average en- 
ergy of the L (R) reservoir. This relation manifests that 
rectification exists if (H^) 7^ {Hr), as in j4^. Specif- 
ically, consider interfaces including 1-dimensional oscil- 
lator chains, = H^"'' + HP°^, where the kinetic en- 
ergy ff^™ is quadratic in momentum, and the potential 
energy per particle is C„g"; n > 2, g is the particle's 
coordinate. In the classical limit using the equipartition 
relation we obtain (i/") = T,j{^ + i). Thermal rectifi- 
cation thus emerges if the reservoirs have a non-identical 
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power n. Note that the separation to three segments (L, 
subsystem, R) is often artificial: The junction could be 
practically made of a single structure with a varying po- 
tential energy, e.g. a nanotube whose composition grad- 
ually changes in space 5l[, leading to inhomogencous 



energy spectra, thus to thermal rectification. 

To conclude, rectification has been obtained here re- 
laying on the DOS asymmetry D^^e) ^ Dfi{e), while 
system- bath interactions are assumed symmetric, \l ~ 
Xn and (7^(0;, e) = gii{Lj,e), energy dependent functions. 
This conclusion is in accord with multitude (numerical) 
observations, demonstrating rectification in two-segment 
dissimilar anharmonic lattices 143. |43|. |44|. |45|. 14611 . 



B. Type-B thermal rectifier 



J{Tl -Tr = ±AT), or by the ratio 7^ = \Jat/J-at\- 
As expected from our general analysis above, a fully har- 
monic system (f^ . a pure spin system (15(11) and a uni- 
form 3LS model (j34p do not rectify heat irrespective of 
the spatial asymmetry embodied in x 7^ 0, since (i) the 
DOS of both contacts are equal, and (ii) subsystem and 
baths are equivalent. This should be emphasized since 
the latter two cases can be viewed as anharmonic struc- 
tures, due to the truncated energy spectra. In contrast, 
for a harmonic bath- TLS-harmonic bath junction (j3ip . 
we obtain the ratio 



7^- 1 



AT 

^ a 



(44) 



Similarly, for a metal- single mode-metal junction, using 
([55)1 at small S, we find 



We explore next the role of the subsystem statistics 
in inducing thermal rectification by further studying the 
inequality ([57)1 . As discussed above, type-A rectifiers 
are constructed assuming that /lCT) ^ fniT), resulting 
from the use of dissimilar reservoirs. However, a more 
careful analysis of (|37p reveals that rectification prevails 
even when /(T) = /^(T) = fniT), as long as Xl ^ Xr 
and the ratio vJ" {—lo) / fL,R{Tu) depends on the temper- 
ature Ty. We show it by rearranging ([57)) . 



n"{-uj) 



1 



1 



^ f{Tc) \Xl 



(43) 



In order for the two sides to deviate, the ratio, e.g., 
{—u))/ /{Th) must depend on the respective temper- 
ature. In other words, the relaxation rates' temperature 
dependence should differ from the central unit particle 
statistics. As shown in Section III.B, the function fv{T) 
reflects the reservoirs statistics, see Eqs. pO)) and (l27l) . 
We therefore classify type-B rectifiers as junctions where 
subsystem and bath differ in their statistics, and the iden- 
tical reservoirs are asymmetrically coupled to the sub- 
system Xl ^ Xn. This observation is in agreement with 
other studies. For example, in Ref. [3^ rectification 
was demonstrated in a quantum-mechanical model where 
photon-mediated heat current flows between two (asym- 
metrically coupled) nonlinear reservoirs. Refs. \M and 
(50| consider a nonlinear subsystem mode while the reser- 
voirs are taken harmonic, yielding rectiflcation due to the 
inclusion of some spatial asymmetry. 

Summarizing, in type-B rectifiers the reservoirs and 
systems-bath couplings are equal at both ends, 0^(6) = 
Df{(e) and gLi^, ^) = 9r{^, ^), but (i) one of the contacts 
is weaker, Xl 7^ Xr, and (ii) the subsystem statistics 
differs from the reservoirs' statistics. 



C. Examples 

In what follows we exemplify type-A and type-B recti- 
fiers. The magnitude of thermal rectification can be es- 
timated by the sum AJ = J at + J~at, where J±at = 



n - 1 



-XfAT. 



(45) 



Thus, in type-B rectifiers the strength of the effect is 
directly linked to the spatial asymmetry, x 7^ [13, HI] . 

We demonstrate next some type-A rectifiers, where the 
reservoirs have distinct properties. In order to simplify 
the presentation we set F = F^ = Tg, i.e. the central 
unit is evenly coupled to the reservoirs. Our first example 
is a phonon bath- HO- spin bath junction, representing 
e.g. an insulating molecule interfacing a dielectric sur- 
face and a nuclear-spin environment. For a schematical 
representation see Fig. [IJa). The current is calculated 
using pT|) with the rates and resulting in 



J = Fw 



n^(cj) -ng(cj) 

, _ ng(-c^) 



The magnitude of thermal rectification is 



nf(-w) l-cxp(-7 



2a) 



-AT' 



n'^{—uj) 1 — exp(- 



2a) 



(46) 



(47) 



-AT I 



It is easy to see that when 



2a; 



3> 1, the rectifier is 



T„±AT 

effectively turned off while for j, "^^at ^ 1 it is operative, 
with 7^ ~ 1 2AT/Ta. 

Our second example is a phonon bath-HO-metal junc- 
tion, representing an electronic to vibrational energy con- 
version device, see Fig. [TJb). This system can be realized 
in a metal-molecule contact, where the vibrational decay 
into the lead (e.g. an STM tip) is bottlenecked by a sin- 
gle vibrational mode. Setting the coupling strength at 
both contacts to be the same, F 



F|^ = F«, we get 



J = ujV- 



1 



1 



5r^ 



(48) 



with the rectification ratio 



n = 



{2 + 5R^){l + 5n^) 



(1 



5r^^ 



1-5r 



Mr ' 



AT 



(49) 
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This expression is similar to (|45p . However, since here 
the reservoirs are distinct (effectively (5^ = 0, J/j ^ 0), 
rectification is obtained for x = 0. 

Fig. m presents a similar instance, vibrations to elec- 
tronic excitations energy exchange, mediated by the ex- 
citation of an anharmonic molecular (TLS) mode. Inter- 
estingly, for small u the solid transfers heat to the metal 
more effectively that the reversed metal-to-dielectric pro- 
cess, reflected by a rectification ratio larger than 1. For 
large spacing the behavior is reversed, and the metal 
better cools down. Consequently, there is an interme- 
diate value (cj ^ 3 for Sp — 0.4) where this nonlinear- 
inhomogeneous junction does not rectify heat. The rec- 
tification ratio obtained here is relatively small, since the 
metal only weakly deviates from the ohmic description 
[39| . This example still demonstrates that by optimizing 
system and interface parameters it may be possible to de- 
sign a junction where the electron bath can be efhciently 
cooled down, while the vibrational energy ineffectively 
dissipates into the metallic bulk, and vice-versa. 

Finally, we examine a spin-subsystem-metal junction. 
In Fig. [5] we show that it can rectify heat in the classical 
limit {lo < T^) while in the quantum regime rectification 
is suppressed. We also modify the system-metal coupling 
strength, and show that it can largely control the rectifi- 
cation ratio (inset). 




0) 



FIG. 4: Electronic to vibrational energy exchange tlirougli a 
TLS mode witli 5r = 0.05 (solid line); Sr = 0.2 (dashed line); 
Sr = 0.4 (dashed-dotted line). Rectification ratio is presented 
as a function of the subsystem energy spacing. Ta = 3, AT — 

1, = 1, r| = r« = 1. 



VI. SUMMARY 

This paper provided a unified description of heat flow 
in two-terminal hybrid structures assuming weak system- 
bath couplings. The underlying origin of nonlinear be- 
havior in various junctions, and its controllability, were 



1 




0.5 1 1.5 2 



CO 

FIG. 5: Spin-HO-metal rectifier (solid line) and a spin-TLS- 
metal rectifier (dashed line). The rectification ratio is pre- 
sented as a function of the subsystem energy spacing. Ta = 
0.5, AT = 0.1, 5r = 0.2, fj,R = 1, Main plot: The subsys- 
tem is equally coupled to the two ends, Vg = Tp = 1. Inset: 
TZ can be tuned by manipulating system-bath interactions, 
= 1, = 0.05. 



explored considering different interfaces: metals, insula- 
tors and noninteracting spins, where the central object 
(subsystem) could represent e.g., a radiation mode or a 
vibrational excitation. 

As a particular example of nonlinear conduction we 
examined thermal rectification. Previous studies of this 
effect were typically based on numerics, utilizing classi- 
cal molecular dynamics tools, where observations were 
deduced within specific molecular force fields 0, IH, 
HIiijHEHl- Here, in contrast, we attempted a gen- 
eral analytical study of the sufficient conditions for the 
onset of thermal rectification in hybrid quantum models. 
We identified two classes of rectifiers: Type-A rectifiers 
where the interfaces are made distinct, and type-B recti- 
fiers where the reservoirs are equivalent, but the system 
and bath quantum statistics differ, in conjunction with 
some spatial (parametric) asymmetry. We note that the 
importance of anharmonicity and asymmetry for mani- 
festing rectification were previously recognized, e.g. in a 
spin-boson model [STj . yet the necessity of the inhomo- 
geneity of the energy spectra was not appreciated. Here 
we clearly observe that a system composed of identical 
anharmonic units cannot rectify heat, unless the energy 
spectra are made non-uniform, in conjimction with some 
spatial asymmetry. 

Our study aims in linking transport characteristics to 
the microscopic Hamiltonian. It might serve as a guide 
for experimentalists pursuing control over energy trans- 
fer in molecular systems tl, 7, 8, 55], and for building 
nanoscale thermal devices [53, |59[ . We have also demon- 
strated that nonlinear thermal effects, and in particu- 
lar thermal rectification, are ubiquitous phenomena that 
could be observed in a variety of systems, phononic [sij . 
electronic [5^] , and photonic [ll|, [s^ . 
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APPENDIX A: Derivation of the heat current 
expression within the master equation formalism 

The aim of this Appendix is to derive the heat current 
expression ([5]) within the Hamiltonian The expecta- 
tion value of the current is formally given by 



J = Tr[Jp]; J - '^[Vl,Hs] + '-[Hs, Vi,]. 
Using the Hamiltonian (|4|) we get 

J — 77 En,mSm.nT^B[^LPn,mBL] 



(Al) 



^ ^ En,mS,n,nT'rB[^RPn,mBR], (A2) 



where = — Tr^ denotes trace over the L 

and R baths, and p is the total density matrix whose 
elements should be calculated in the long time limit, since 
expression (jAl[) is valid in steady-state situations only 
[53j . We find it useful to rearrange Eq. (|A2p as follows. 



J — T: -E'n^mS'm^nTrB [Ai(p„ 



nyra 

(A3) 

We would like to express the current in terms of the long- 
time population P„ . The Liouville equation of motion for 



Pn.m 



rnP' 



n.mh'n.m 



n,pPp,m 
u p 

(A4) 

with the transformed operators B,y{t) = e^^"^ B,je~^^"^; 
V = L,R. The index p counts the subsystem states. This 
equation can be formally integrated to yield 

(r)S,(T)]dT. (A5) 

The initial condition was already neglected since it will 
not contribute after tracing out the bath in (|A3|) . as- 
suming TrB[i3(i), p(0)] = 0, i.e. the mean value of 



the interaction Hamiltonian, averaged over the initial 
density matrix, is zero. Plugging (|A5[) into (|A3p we 
obtain order terms. Factorizing the density ma- 
trix at all times, p{t) = pL{TL)pR{TR)(T{t), py{Ty) = 
e~""/T^ /Ti,[e-"°/^"]- a{t) = Trs[p(i)], and taking the 
markovian limit, replacing (t(t) — > a{t), we get 



TrB[pn,m{t)BL(t)] 



[BL{t)BL{T))T, ^p,m 



{BL{T)BL{t))TL(7n,p{t)Sp,m cLt. 



(A6) 



Here (S^(r)B^(0))T„ = Tr^[p^(r^)B^(r)B^(0)] is the 
correlation function of the v environment. In de- 
riving (jA6[) we disregarded mixed terms of the form 
{BL(t)Bfi{T)), since the terminals are not correlated. 
Our next assumption is that coherences are negligi- 
ble in the weak-coupling scheme, given the preparation 
cTn=^m{0) ~ 0, thus wc disregard all the nondiagonal terms 
in (jA6p . Further, performing the transformation x — r—t 
and extending the integral to infinity (markovian approx- 
imation), we obtain 



TrB[Pn,m{t)BL{t)] 
1-0 



-iX 



L I e 

— OO 



[Bl{0)Bl{x))t,, {t)Sn 



-(BLix)BLiO))TLa„^nit)Sn,m dx 

Similarly, we find that 

TrB[p,„,„(t)Si(i)] « 



(A7) 



lo 



{BL{x)BLmT, {t)S 

n 



{Bl{0)Bl{x))t, .n 



dx. 



(A8) 



Combining the last two expressions, we come by (using 
Tr_B[(p„^„(i) - Pm,7i{t))BL{t)] w 



iXrSn 



Pn{t) I e'''--%BL{x)BLmT, 



'{BL{0)BL{x)h 



(A9) 



We return now to the heat current expression (|A3p . and 
identify the v-hath induced relaxation rates by 



dre 



B,{t)B,{Q))j, .(AlO) 



Making use of (|A9|) . conjoined with the analogous R- 
bath expression, we get a weak-coupling expression for 
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the heat current, defined positive when flowing L to i?, 

J — ~2 En^m\Sm,n\'^\Pnkn^m{T^L) ~ -Pmfem^n )] 

+ ^ E En,rn\SrnA''[Pnk^-.rn{TR) - Prnk^^^iTn)], 
n'>rn 

(All) 

or more compactly 

(A12) 



This is our final result, a second-order (A^) expression for 
the heat current, obtained within the master equation 
approach. The populations P„ should be calculated at 
long time, using ©. 
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